Data are loaded from .Rdata.
setwd(getwd())
load("tms_data_preprocessed.Rdata")
bname = "tms_analyses"
models_task <- list(
formula(probe.response ~ zbv * zlog.apen + (1|subj/condition)),
formula(probe.response ~ zbv * zlog.apen + probeix + (1|subj/condition)),
formula(probe.response ~ zbv * zlog.apen + probeix + block_num + (1|subj/condition)),
formula(probe.response ~ zbv * zlog.apen + probeix + block_num + condition + (1|subj/condition)),
formula(probe.response ~ zbv * zlog.apen + probeix + block_num + condition + randomization + (1|subj/condition))
)
descriptions_task=c(
"BV x AE",
"BV x AE + trial",
"BV x AE + trial + block",
"BV x AE + trial + block + condition",
"BV x AE + trial + block + condition + randomization"
)
loos_task <- load.cache.var("loos_task", base = bname) # load loos
CACHE> loading loos_task from cache/vars/tms_analyses_loos_task.RData
mod_selection_res <- as.data.frame(loos_task$diffs) %>% select(c(elpd_diff, se_diff)) # from loos select elpd_diff and se_diff, convert to data frame
models_elpd_diff <- c(
"BV x AE + probe + block",
"BV x AE + probe + block + condition",
"BV x AE + probe + block + condition + randomization",
"BV x AE + probe",
"BV x AE"
) #List models in the same order as in elpd_diff. Those are labels for the graphical table
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff) # create a data frame with model labels and loo diagnostics
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff") # name columns appropriately
mod_selection_res = mod_selection_res %>%
mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>%
rename_all(~gsub("\\.", " ", .)) # format table
alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l")) # align rows
mod_selection_res %>% # plot model selection table
kbl(caption="LOO-CV: MW", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(2, bold = TRUE, color = "red")
| Model | elpd_diff | se_diff |
|---|---|---|
| BV x AE + probe + block | 0.00 | 0.00 |
| BV x AE + probe + block + condition | -0.24 | 1.25 |
| BV x AE + probe + block + condition + randomization | -0.29 | 1.26 |
| BV x AE + probe | -13.58 | 5.44 |
| BV x AE | -17.27 | 6.33 |
color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))
mod_task03 <- load.cache.var("mod_task03", base = bname)
CACHE> loading mod_task03 from cache/vars/tms_analyses_mod_task03.RData
mcmc_intervals(mod_task03, prob_outer = 0.95, pars = c("b_zbv",
"b_zlog.apen",
"b_block_num",
"b_probeix",
"b_zbv:zlog.apen",
# "b_visit2",
"b_conditionsham_rhTMS",
"b_conditionsham_arrhTMS",
"b_conditionactive_rhTMS",
"b_conditionactive_arrhTMS"
)) +
ggplot2::scale_y_discrete(labels = c("b_zbv" = "BV",
"b_zlog.apen" = "AE",
"b_block_num" = "Block",
"b_probeix" = "Probe number",
"b_zbv:zlog.apen" = "BV x AE",
# "b_visit2" = "Visit",
"b_conditionsham_rhTMS" = "Sham rhTMS",
"b_conditionsham_arrhTMS" = "Sham arrhTMS",
"b_conditionactive_rhTMS" = "Active rhTMS",
"b_conditionactive_arrhTMS" = "Active arrhTMS")) +
geom_vline(xintercept = 0,color = "red", size=0.75) + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black','black', 'black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: MW")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Fig. 1. Model 7: On-task Score: coefficients for all parameters of interest.
mw_hyp = hypothesis(mod_task03, c("zbv<0", "zlog.apen>0", "block_num<0", "probeix<0", "zbv:zlog.apen<0", "conditionactive_rhTMS>0", "conditionsham_rhTMS>0", "conditionactive_arrhTMS>0", "conditionsham_arrhTMS>0"), alpha = 0.025) # test hypotheses for model 03
mw_hyp_table <- mw_hyp$hypothesis %>% select(-Star) # select all columns from the resulting object except for Star
alignment= map_chr(mw_hyp_table, ~ifelse(class(.x)=="numeric", "r","l")) # alogn rows
hypotheses <- c("BV < 0", "AE > 0", "Block < 0", "Probe number < 0", "BV x AE< 0", "Active rhTMS > 0", "Sham rhTMS > 0", "Active arrhTMS > 0", "Sham arrhTMS > 0") # labels for the hypothesis column
mw_hyp_table = mw_hyp_table %>% # table formatting
mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>%
rename_all(~gsub("\\.", " ", .))
mw_hyp_table %>% # plot table
kbl(caption="Results: MW", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(7, 8, 9), bold = TRUE) %>%
row_spec(6, bold = TRUE, color = "red")
| Hypothesis | Estimate | Est Error | CI Lower | CI Upper | Evid Ratio | Post Prob |
|---|---|---|---|---|---|---|
| BV < 0 | -0.03 | 0.05 | -0.13 | 0.06 | 3.37 | 0.77 |
| AE > 0 | 0.06 | 0.05 | -0.03 | 0.15 | 9.58 | 0.91 |
| Block < 0 | -0.11 | 0.02 | -0.15 | -0.07 | Inf | 1.00 |
| Probe number < 0 | -0.04 | 0.01 | -0.07 | -0.02 | 713.29 | 1.00 |
| BV x AE< 0 | -0.04 | 0.05 | -0.14 | 0.06 | 3.75 | 0.79 |
| Active rhTMS > 0 | 0.28 | 0.26 | -0.23 | 0.79 | 6.52 | 0.87 |
| Sham rhTMS > 0 | -0.23 | 0.26 | -0.75 | 0.26 | 0.22 | 0.18 |
| Active arrhTMS > 0 | 0.06 | 0.26 | -0.44 | 0.58 | 1.43 | 0.59 |
| Sham arrhTMS > 0 | -0.32 | 0.26 | -0.84 | 0.19 | 0.12 | 0.11 |
tms_data.nback_z <- tms_data.nback %>%
reshape2::melt(id.vars = c( "focus", "probe.response", "intention", "somnolence", "condition", "visit", "subj"), measure.vars = c("zlog.apen", "zbv"), varnames = c("Variable", "Score"))
data.table::setnames(tms_data.nback_z, old = c("focus", "probe.response", "intention", "somnolence", "condition", "visit", "subj",'variable', "value"), new = c('Focus', "On-task Score", "Intention", "Alertness", "Condition", "Visit", "Subject",'Measure', "Z-score"))
Fig. 2. MW Mean ± 1SE across conditions.
loos_bv=load.cache.var("loos_bv", base=bname)
CACHE> loading loos_bv from cache/vars/tms_analyses_loos_bv.RData
# same procedure for BV and AE
as.data.frame(loos_bv$diffs)
mod_selection_res <- as.data.frame(loos_bv$diffs) %>% select(c(elpd_diff, se_diff))
models_elpd_diff <- c(
"probe + block + condition + visit",
"probe + block",
"probe + block + condition",
"probe"
)
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff)
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff")
mod_selection_res = mod_selection_res %>%
mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>%
rename_all(~gsub("\\.", " ", .))
alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l"))
elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic
mod_bv03 0.000000 0.000000 -749.5354 37.47914 44.65426 1.245530 1499.071
mod_bv01 -2.226843 4.354602 -751.7623 37.03687 40.26908 1.058170 1503.525
mod_bv02 -3.282269 3.992372 -752.8177 36.94252 41.68978 1.117909 1505.635
mod_bv00 -4.509851 5.198430 -754.0453 37.10229 38.84448 1.067001 1508.091
se_looic
mod_bv03 74.95829
mod_bv01 74.07375
mod_bv02 73.88504
mod_bv00 74.20458
mod_selection_res %>%
kbl(caption="LOO-CV: BV", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(1, bold = TRUE, color = "red")
| Model | elpd_diff | se_diff |
|---|---|---|
| probe + block + condition + visit | 0.00 | 0.00 |
| probe + block | -2.23 | 4.35 |
| probe + block + condition | -3.28 | 3.99 |
| probe | -4.51 | 5.20 |
color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))
mod_bv03 <- load.cache.var("mod_bv03", base = bname)
mcmc_intervals(mod_bv03, prob_outer = 0.95, pars = c(
"b_block_num",
"b_probeix",
"b_visit2",
"b_conditionsham_rhTMS",
"b_conditionsham_arrhTMS",
"b_conditionactive_rhTMS",
"b_conditionactive_arrhTMS"
)) +
ggplot2::scale_y_discrete(labels = c("b_block_num" = "Block",
"b_probeix" = "Probe number",
"b_visit2" = "Visit",
"b_conditionsham_rhTMS" = "Sham rhTMS",
"b_conditionsham_arrhTMS" = "Sham arrhTMS",
"b_conditionactive_rhTMS" = "Active rhTMS",
"b_conditionactive_arrhTMS" = "Active arrhTMS")) +
geom_vline(xintercept = 0,color = "red", size=0.75) + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: BV")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Fig. 6. BV: coefficients for all parameters of interest.
CACHE> loading mod_bv03 from cache/vars/tms_analyses_mod_bv03.RData
bv_hyp = hypothesis(mod_bv03, c("block_num>0", "probeix>0", "conditionactive_rhTMS>0", "conditionsham_rhTMS>0", "conditionactive_arrhTMS>0", "conditionsham_arrhTMS>0"), alpha = 0.025)
bv_hyp_table <- bv_hyp$hypothesis %>% select(-Star)
alignment= map_chr(bv_hyp_table, ~ifelse(class(.x)=="numeric", "r","l"))
hypotheses <- c("Block > 0", "Probe number > 0", "Active rhTMS > 0", "Sham rhTMS > 0", "Active arrhTMS > 0", "Sham arrhTMS > 0")
bv_hyp_table = bv_hyp_table %>%
mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>%
rename_all(~gsub("\\.", " ", .))
bv_hyp_table %>%
kbl(caption="Results: BV", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(4, 5, 6), bold = TRUE) %>%
row_spec(3, bold = TRUE, color = "red")
| Hypothesis | Estimate | Est Error | CI Lower | CI Upper | Evid Ratio | Post Prob |
|---|---|---|---|---|---|---|
| Block > 0 | -0.02 | 0.01 | -0.04 | 0.00 | 0.01 | 0.01 |
| Probe number > 0 | 0.00 | 0.01 | -0.01 | 0.01 | 1.17 | 0.54 |
| Active rhTMS > 0 | 0.18 | 0.09 | 0.01 | 0.35 | 50.02 | 0.98 |
| Sham rhTMS > 0 | 0.04 | 0.09 | -0.12 | 0.22 | 2.36 | 0.70 |
| Active arrhTMS > 0 | 0.13 | 0.09 | -0.04 | 0.30 | 13.90 | 0.93 |
| Sham arrhTMS > 0 | 0.04 | 0.08 | -0.12 | 0.21 | 2.21 | 0.69 |
loos_apen=load.cache.var("loos_apen", base=bname)
CACHE> loading loos_apen from cache/vars/tms_analyses_loos_apen.RData
mod_selection_res <- as.data.frame(loos_apen$diffs) %>% select(c(elpd_diff, se_diff))
models_elpd_diff <- c(
"probe + block",
"probe",
"probe + block + condition",
"probe + block + condition + visit"
)
mod_selection_res <- data.frame(models_elpd_diff, mod_selection_res$elpd_diff, mod_selection_res$se_diff)
names(mod_selection_res) <- c("Model", "elpd_diff", "se_diff")
mod_selection_res = mod_selection_res %>%
mutate(across(where(is.numeric), ~comma(., accuracy=0.01))) %>%
rename_all(~gsub("\\.", " ", .))
alignment= map_chr(mod_selection_res, ~ifelse(class(.x)=="numeric", "r","l"))
mod_selection_res %>%
kbl(caption="LOO-CV: AE", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(4, bold = TRUE, color = "red")
| Model | elpd_diff | se_diff |
|---|---|---|
| probe + block | 0.00 | 0.00 |
| probe | -0.77 | 1.99 |
| probe + block + condition | -0.97 | 1.97 |
| probe + block + condition + visit | -1.45 | 2.01 |
color_scheme_set("viridisE")
bayesplot_theme_set(theme_bw(base_size = 15, base_family = "sans"))
mod_apen03 <- load.cache.var("mod_apen03", base = bname)
mcmc_intervals(mod_apen03, prob_outer = 0.95, pars = c(
"b_block_num",
"b_probeix",
"b_visit2",
"b_conditionsham_rhTMS",
"b_conditionsham_arrhTMS",
"b_conditionactive_rhTMS",
"b_conditionactive_arrhTMS"
)) +
ggplot2::scale_y_discrete(labels = c("b_block_num" = "Block",
"b_probeix" = "Probe number",
"b_visit2" = "Visit",
"b_conditionsham_rhTMS" = "Sham rhTMS",
"b_conditionsham_arrhTMS" = "Sham arrhTMS",
"b_conditionactive_rhTMS" = "Active rhTMS",
"b_conditionactive_arrhTMS" = "Active arrhTMS")) +
geom_vline(xintercept = 0,color = "red", size=0.75) + theme(axis.text.y = element_text(size=18, face = "plain", colour = c('black', 'black', 'black', 'black', 'black', 'red', 'black'))) + ggtitle("Model Coefficients: AE")
Scale for 'y' is already present. Adding another scale for 'y', which will
replace the existing scale.
Warning: Vectorized input to `element_text()` is not officially supported.
Results may be unexpected or may change in future versions of ggplot2.
CACHE> loading mod_apen03 from cache/vars/tms_analyses_mod_apen03.RData
Fig. 7. AE: coefficients for all parameters of interest.
apen_hyp = hypothesis(mod_apen03, c("block_num<0", "probeix<0", "conditionactive_rhTMS<0", "conditionsham_rhTMS<0", "conditionactive_arrhTMS<0", "conditionsham_arrhTMS<0"), alpha = 0.025)
apen_hyp_table <- apen_hyp$hypothesis %>% select(-Star)
alignment= map_chr(apen_hyp_table, ~ifelse(class(.x)=="numeric", "r","l"))
hypotheses <- c("Block < 0", "Probe number < 0", "Active rhTMS < 0", "Sham rhTMS < 0", "Active arrhTMS < 0", "Sham arrhTMS < 0")
apen_hyp_table = apen_hyp_table %>%
mutate(across(where(is.numeric), ~comma(., accuracy=0.01)), Hypothesis = hypotheses) %>%
rename_all(~gsub("\\.", " ", .))
apen_hyp_table %>%
kbl(caption="Results: AE", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(4, 5, 6), bold = TRUE) %>%
row_spec(3, bold = TRUE, color = "red")
| Hypothesis | Estimate | Est Error | CI Lower | CI Upper | Evid Ratio | Post Prob |
|---|---|---|---|---|---|---|
| Block < 0 | -0.02 | 0.01 | -0.05 | 0.00 | 23.81 | 0.96 |
| Probe number < 0 | 0.02 | 0.01 | 0.00 | 0.04 | 0.01 | 0.01 |
| Active rhTMS < 0 | -0.27 | 0.12 | -0.50 | -0.04 | 83.03 | 0.99 |
| Sham rhTMS < 0 | -0.22 | 0.12 | -0.45 | 0.01 | 32.78 | 0.97 |
| Active arrhTMS < 0 | -0.23 | 0.12 | -0.47 | -0.01 | 43.25 | 0.98 |
| Sham arrhTMS < 0 | -0.12 | 0.12 | -0.35 | 0.11 | 5.97 | 0.86 |
# Check data normality assumption. Spoiler alert: violated.
shapiro.test(tms_data.nback$probe.response)
shapiro.test(tms_data.nback$zlog.apen)
shapiro.test(tms_data.nback$bv)
Shapiro-Wilk normality test
data: tms_data.nback$probe.response
W = 0.87298, p-value < 2.2e-16
Shapiro-Wilk normality test
data: tms_data.nback$zlog.apen
W = 0.9406, p-value < 2.2e-16
Shapiro-Wilk normality test
data: tms_data.nback$bv
W = 0.5741, p-value < 2.2e-16
Distributions are non-normal.
kruskal.test(probe.response ~ condition, data = tms_data.nback) # groups are different
Kruskal-Wallis rank sum test
data: probe.response by condition
Kruskal-Wallis chi-squared = 15.76, df = 4, p-value = 0.003359
tms_data.nback$condition <- fct_relevel(tms_data.nback$condition, "baseline", "active_rhTMS", "active_arrhTMS", "sham_rhTMS", "sham_arrhTMS") # relevel condition factors
tms_comparison_list <- list(c("baseline", "active_rhTMS"), c("baseline", "active_arrhTMS"), c("baseline", "sham_rhTMS"), c("baseline", "sham_arrhTMS"), c("active_rhTMS", "active_arrhTMS"), c("active_rhTMS", "sham_rhTMS"), c("active_rhTMS", "sham_arrhTMS"), c("active_arrhTMS", "sham_rhTMS"), c("active_arrhTMS", "sham_arrhTMS"), c("sham_rhTMS", "sham_arrhTMS")) # define the exact comparison pairs
pairwise_wilcox_results <- data.frame(matrix(ncol = 4, nrow = 10)) # create an empty data frame for Wilcoxon results
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value", "Effect size, r") # name the columns
pairwise_wilcox_results$Comparison <- tms_comparison_list # assign comparisons to the corresponding column
eff_sizes_mw <- rstatix::wilcox_effsize(tms_data.nback, probe.response ~ condition, paired = FALSE) # compute effect size estimates
i = 1
for (comparison in tms_comparison_list) { # perform Wilcox test for each comparison pair
r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$probe.response, filter(tms_data.nback, condition == comparison[2])$probe.response, paired = FALSE)
pairwise_wilcox_results[i, 2:4] <- rbind(bind_cols(r$statistic, r$p.value, eff_sizes_mw$effsize[[i]])) #bind rows
i = i + 1
}
comparisons <- c("baseline vs. active rhTMS", "baseline vs. active arrhTMS", "baseline vs. sham rhTMS", "baseline vs. sham arrhTMS", "active rhTMS vs. active arrhTMS", "active rhTMS vs. sham rhTMS", "active rhTMS vs. sham arrhTMS", "active arrhTMS vs. sham rhTMS", "active arrhTMS vs. sham arrhTMS", "sham rhTMS vs. sham arrhTMS") # labels for comparisons
pairwise_wilcox_results <- pairwise_wilcox_results %>% # format table
mutate(across("p-value", ~comma(., accuracy=0.01)), across("Effect size, r", ~comma(., accuracy=0.01)),Comparison = comparisons) %>%
rename_all(~gsub("\\.", " ", .))
pairwise_wilcox_results %>% # plot table
kbl(caption="Pairwise Wilcoxon test: MW", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(3, 4, 6, 7), bold = TRUE) %>%
row_spec(1, color = "red")
| Comparison | W-statistic | p-value | Effect size, r |
|---|---|---|---|
| baseline vs. active rhTMS | 18312 | 0.41 | 0.04 |
| baseline vs. active arrhTMS | 20483 | 0.23 | 0.06 |
| baseline vs. sham rhTMS | 21615 | 0.02 | 0.11 |
| baseline vs. sham arrhTMS | 22196 | 0.00 | 0.14 |
| active rhTMS vs. active arrhTMS | 14191 | 0.07 | 0.10 |
| active rhTMS vs. sham rhTMS | 14899 | 0.01 | 0.15 |
| active rhTMS vs. sham arrhTMS | 15308 | 0.00 | 0.18 |
| active arrhTMS vs. sham rhTMS | 13626 | 0.29 | 0.06 |
| active arrhTMS vs. sham arrhTMS | 13944 | 0.14 | 0.08 |
| sham rhTMS vs. sham arrhTMS | 13088 | 0.71 | 0.02 |
# do the same for BV and AE
kruskal.test(zlog.apen ~ condition, data = tms_data.nback) # groups are different
Kruskal-Wallis rank sum test
data: zlog.apen by condition
Kruskal-Wallis chi-squared = 9.6253, df = 4, p-value = 0.04724
pairwise_wilcox_results <- data.frame(matrix(ncol = 4, nrow = 10))
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value", "Effect size, r")
pairwise_wilcox_results$Comparison <- tms_comparison_list
eff_sizes_mw <- rstatix::wilcox_effsize(tms_data.nback, zlog.apen ~ condition, paired = FALSE)
i = 1
for (comparison in tms_comparison_list) {
r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$zlog.apen, filter(tms_data.nback, condition == comparison[2])$zlog.apen, paired = FALSE)
pairwise_wilcox_results[i, 2:4] <- rbind(bind_cols(r$statistic, r$p.value, eff_sizes_mw$effsize[[i]]))
i = i + 1
}
pairwise_wilcox_results <- pairwise_wilcox_results %>%
mutate(across("p-value", ~comma(., accuracy=0.01)), across("Effect size, r", ~comma(., accuracy=0.01)),Comparison = comparisons) %>%
rename_all(~gsub("\\.", " ", .))
pairwise_wilcox_results %>%
kbl(caption="Pairwise Wilcoxon test: AE", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(1, 2), bold = TRUE) %>%
row_spec(1, color = "red")
| Comparison | W-statistic | p-value | Effect size, r |
|---|---|---|---|
| baseline vs. active rhTMS | 22138.5 | 0.01 | 0.13 |
| baseline vs. active arrhTMS | 21689.0 | 0.03 | 0.11 |
| baseline vs. sham rhTMS | 21311.5 | 0.06 | 0.09 |
| baseline vs. sham arrhTMS | 19889.5 | 0.54 | 0.03 |
| active rhTMS vs. active arrhTMS | 12561.5 | 0.77 | 0.02 |
| active rhTMS vs. sham rhTMS | 12431.0 | 0.66 | 0.02 |
| active rhTMS vs. sham arrhTMS | 11337.5 | 0.08 | 0.10 |
| active arrhTMS vs. sham rhTMS | 12673.5 | 0.88 | 0.01 |
| active arrhTMS vs. sham arrhTMS | 11567.5 | 0.14 | 0.08 |
| sham rhTMS vs. sham arrhTMS | 11767.5 | 0.21 | 0.07 |
pairwise_wilcox_results <- data.frame(matrix(ncol = 3, nrow = 10))
colnames(pairwise_wilcox_results) <- c("Comparison", "W-statistic", "p-value")
pairwise_wilcox_results$Comparison <- tms_comparison_list
i = 1
for (comparison in tms_comparison_list) {
r <- wilcox.test(filter(tms_data.nback, condition == comparison[1])$zbv, filter(tms_data.nback, condition == comparison[2])$zbv, paired = FALSE)
pairwise_wilcox_results[i, 2:3] <- rbind( bind_cols(r$statistic, r$p.value))
i = i + 1
}
pairwise_wilcox_results <- pairwise_wilcox_results %>%
mutate(across("p-value", ~comma(., accuracy=0.01)), Comparison = comparisons) %>%
rename_all(~gsub("\\.", " ", .))
pairwise_wilcox_results %>%
kbl(caption="Pairwise Wilcoxon test: BV", align=alignment) %>%
kable_classic(full_width=TRUE, html_font="Times") %>%
row_spec(c(2), bold = TRUE) %>%
row_spec(1, color = "red")
| Comparison | W-statistic | p-value |
|---|---|---|
| baseline vs. active rhTMS | 16467 | 0.02 |
| baseline vs. active arrhTMS | 16813 | 0.04 |
| baseline vs. sham rhTMS | 18209 | 0.38 |
| baseline vs. sham arrhTMS | 17964 | 0.28 |
| active rhTMS vs. active arrhTMS | 13154 | 0.67 |
| active rhTMS vs. sham rhTMS | 13859 | 0.20 |
| active rhTMS vs. sham arrhTMS | 13786 | 0.23 |
| active arrhTMS vs. sham rhTMS | 13625 | 0.32 |
| active arrhTMS vs. sham arrhTMS | 13615 | 0.32 |
| sham rhTMS vs. sham arrhTMS | 12678 | 0.88 |
wilcox.test(filter(tms_data.nback_z, (`Condition` == "baseline" & `Measure` == "zlog.apen"))$`Z-score`,
filter(tms_data.nback_z, (`Condition` == "active_rhTMS" & `Measure` == "zlog.apen"))$`Z-score`,
paired = FALSE)
rstatix::wilcox_effsize(tms_data.nback, zlog.apen ~ condition, paired = FALSE)
Wilcoxon rank sum test with continuity correction
data: filter(tms_data.nback_z, (Condition == "baseline" & Measure == "zlog.apen"))$`Z-score` and filter(tms_data.nback_z, (Condition == "active_rhTMS" & Measure == "zlog.apen"))$`Z-score`
W = 22138, p-value = 0.009497
alternative hypothesis: true location shift is not equal to 0
# A tibble: 10 x 7
.y. group1 group2 effsize n1 n2 magnitude
* <chr> <chr> <chr> <dbl> <int> <int> <ord>
1 zlog.apen baseline active_rhTMS 0.130 240 160 small
2 zlog.apen baseline active_arrhTMS 0.110 240 160 small
3 zlog.apen baseline sham_rhTMS 0.0932 240 160 small
4 zlog.apen baseline sham_arrhTMS 0.0304 240 160 small
5 zlog.apen active_rhTMS active_arrhTMS 0.0161 160 160 small
6 zlog.apen active_rhTMS sham_rhTMS 0.0249 160 160 small
7 zlog.apen active_rhTMS sham_arrhTMS 0.0988 160 160 small
8 zlog.apen active_arrhTMS sham_rhTMS 0.00855 160 160 small
9 zlog.apen active_arrhTMS sham_arrhTMS 0.0833 160 160 small
10 zlog.apen sham_rhTMS sham_arrhTMS 0.0697 160 160 small
kruskal.test(zbv ~ condition, data = tms_data.nback) # groups are NOT different
Kruskal-Wallis rank sum test
data: zbv by condition
Kruskal-Wallis chi-squared = 7.4528, df = 4, p-value = 0.1138
tms_data.nback_z %>%
ggplot(aes(x= `Condition`, y =`Z-score`, group = `Measure`, shape = `Condition`, color = `Measure`)) +
geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
theme_set(theme_bw()) +
theme(text = element_text(size = 15)) +
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=0.4))
Fig. 3. Mean ± 1SE changes of AE and BV across conditions.
tms_data.nback_z %>%
ggplot(aes(x= `On-task Score`, y =`Z-score`, group = `Measure`, color = `Measure`)) +
geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
facet_wrap(~`Condition`) +
scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
theme(text = element_text(size = 12))
Fig. 4. AE & BV versus On-task score: Mean ± 1SE across conditions. Our subjects are either bad at tapping with the metronome or misunderstood this part of the task.
tms_data.nback_z %>%
ggplot(aes(x= `Focus`, y =`Z-score`, group = `Measure`, color=`Measure`)) +
geom_pointrange(aes(shape = `Condition`),stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
facet_wrap(~ `Condition`)
Fig. 5. AE & BV versus On-task score: BV x AE interaction versus on-off task states. acrive rhTMS seems to flip the interaction: subjects have worse performance when they think they are on task and vice-versa.
tms_data.nback %>%
ggplot(aes(x = zbv)) +
geom_histogram(binwidth = 0.1) +
labs(y = "Count", x = "Z-score: BV") +
theme(text = element_text(size = 18))
tms_data.nback %>%
ggplot(aes(x = zlog.apen)) +
geom_histogram(binwidth = 0.1) +
labs(y = "Count", x = "Z-score: AE") +
theme(text = element_text(size = 18))
tms_data.nback %>%
ggplot(aes(x = probe.response)) +
geom_histogram(binwidth = 0.5) +
labs(y = "Count", x = "On-task score") +
theme(text = element_text(size = 18))
AE and BV distributions
tms_data.nback_z %>%
ggplot(aes(x= `Condition`, y =`On-task Score`, group = 1)) +
geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
facet_wrap(~ `Subject`, scales = "free", ncol = 2, as.table = FALSE) +
theme_set(theme_bw()) +
geom_vline(xintercept = "active_rhTMS", color = "red", size=0.5, linetype= "dotdash") +
theme(axis.text.x = element_text(angle = 20, vjust = 0.5, hjust=0.4)) +
theme(plot.margin = margin(1,1,1.5,1.2, "cm"))
AE and BV distributions
tms_data.nback_z %>%
ggplot(aes(x= `On-task Score`, y =`Z-score`, group = `Measure`, color = `Measure`)) +
geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
facet_wrap(~`Condition`) +
scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
theme(text = element_text(size = 15))
tms_data.nback_z %>%
ggplot(aes(x= `Focus`, y =`Z-score`, group = `Measure`, color=`Measure`)) +
geom_pointrange(stat="summary", fun.data=mean_se, fun.args = list(mult=1), position=position_dodge(0.05)) +
geom_line(stat="summary", fun.data=mean_se, fun.args = list(mult=2)) +
scale_color_manual(labels = c("AE", "BV"), values = c("blue", "red")) +
facet_wrap(~ `Condition`)+
theme(text = element_text(size = 18))